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Although Laplace's equation is simple, the region over which it is 
to be solved is often complicated. Both the shape of the region and the 
boundary conditions can induce solutions <t> which are singular at iso- 
lated points on the boundary of the region. 

Boundary integral equation methods are well-suited to the problem, 
reducing a two-dimensional partial differential equation to a one- 
dimensional integral equation. Unfortunately, the standard boundary 
integral equation methods lead to an ill-conditioned set of linear 
equations, restricting the achievable accuracy in the approximate so- 
lution. 

This paper describes an improved boundary integral method. A new 
integral equation is derived. Laplace's equation is reduced to solving 
two coupled, one-dimensional integral equations. The resulting linear 
equations are well-conditioned. 

A program package for solving Laplace's equation has been devel- 
oped. The package solves Laplace's equation in two dimensions or in 
three dimensions with axial symmetry. The region may extend to in- 
finity, and may be multiply-connected. In addition to smooth basis 
functions, the program automatically includes appropriate singular 
basis functions, greatly improving the achievable accuracy for regions 
with corners. 



I. INTRODUCTION 

Laplace's equation frequently arises in modeling physical problems, 
especially in electromagnetism, in thermal flow, and in fluid flow. In two 
dimensions, Laplace's equation is 

d 2 $ d 2 $ _ 

a* 2 dy* ' 

2797 




Fig. 1 — Region used in an analysis of an electrostatic lens. 

and in three dimensions, 

d 2 <i> d 2 3> d 2 «i> 
dx 2 dy 2 dz 2 

To complete the specification of a particular problem, a region on which 
to solve Laplace's equation must be specified, plus boundary conditions 
on the boundary of the region. 

As compensation for the simplicity of the partial differential equation, 
the region over which Laplace's equation is to be solved is often com- 
plicated. Figure 1 shows the region used by the author in an unpublished 
analysis of an electrostatic lens. The solution is singular at the re-entrant 
corners. (By singular, we mean that $ has a finite limit as the corner is 
approached, but that some derivatives of $ do not have a finite limit.) 
The singularity is a consequence of the region itself, not of any particular 
boundary conditions. In fact, the solution is singular unless very special 
boundary conditions are prescribed. 

Even with a rectangular region, the solution can be singular at isolated 
points. Figure 2 is an example, a thin-film capacitor with metal top and 
bottom contacts. To obtain its capacitance, Laplace's equation must be 
solved inside the rectangle. The boundary conditions are $ = 1 on the 
top contact, $ = on the bottom contact, and zero normal derivative, 
d$/dn = 0, on the remainder of the boundary. (The definition of the 
normal derivative is given in the next section.) At the center edge of the 
top contact, the solution is singular. 

Standard methods for elliptic partial differential equations include 
finite difference and finite element methods. Both methods require a 
grid, usually rectangular or triangular, everywhere inside the region. 
Thus the region must be bounded. Both methods are difficult to apply 
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Fig. 2 — Region used in an analysis of a thin-film capacitor. The potential is singular 
at the center of the top side. 

to complicated regions; if the true solution has singularities, accuracy 
is usually poor unless heroic measures are taken. Neither method is 
suitable for a package for general regions and boundary conditions. 

Laplace's equation is the simplest elliptic partial differential equation, 
and has been the subject of a great deal of analysis. Special methods for 
Laplace's equation are available, methods that do not work for general 
elliptic partial differential equations. 

Special methods for Laplace's equation include the so-called "fast 
Poisson solvers." * They can quickly solve V 2 <j>(x,y) = f{x,y) if the region 
and boundary conditions are sufficiently simple. However, even Fig. 2 
is not simple enough because of the mixed boundary conditions on the 
top boundary. The fast Poisson solvers have great utility for special 
problems, but are not appropriate for a general Laplace package. Recent 
research (Ref. 2, for example) indicates how these methods may be ex- 
tended in the future. 

1. 1 The boundary Integral equation method 

The most useful special method for Laplace's equation is the boundary 
integral equation method. The basic method has been known for many 
years, 3,4 but has enjoyed a renewed popularity since the advent of large 
digital computers. A few representative references are Refs. 5 to 9. A 
two-dimensional partial differential equation is reduced to a one-di- 
mensional integral equation. Similarly, a three-dimensional partial 
differential equation can be reduced to a two-dimensional integral 
equation. The integral equation involves only the geometry and the 
values of $ and d$/dn on the boundary. Multiply-connected regions pose 
no added difficulty. After the integral equation has been solved ap- 
proximately, another integral can be done to evaluate V<i> and <I> at any 
point inside the region. 

The boundary integral equation method has been quite successful, 
providing fast and inexpensive solutions for Laplace's equation in two 
dimensions. The usual implementation does have several difficulties. 
First, the integral equation is a Fredholm integral equation of the first 
kind for d<J>/dn, and consequently is ill-conditioned. (For either a Diri- 
chlet or a Neumann problem, a well-conditioned integral equation is 
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Fig. 3 — Illustration of the definitions of s, x s , y s , n„, and F. 

available, but not for mixed boundary conditions.) Matrices generated 
for approximate solutions to the integral equation are ill-conditioned, 
and d$/dn cannot be found accurately. Second, three-dimensional 
problems with axial symmetry are essentially two-dimensional problems, 
but no provision is made for their solution. (The next two objections do 
not apply to Ref. 9.) Third, the unknown $ and d$/dn are approximated 
by low-order polynomials on sections of the boundary. Convergence 
requires many coefficients if accuracy of more than a few percent is re- 
quired. Finally, no provision is provided for dealing with singu- 
larities. 

The present paper describes an improved boundary integral method 
which counters all the above difficulties. Two coupled integral equations 
are used; the combination leads to a well-conditioned matrix. Both $ 
and d$/dn can be obtained accurately. Higher-order approximations 
for the unknown $ and d$/dn are used. Corner singularities are recog- 
nized automatically, and special approximating functions are used. 
Axisymmetric problems are solved by the same program. Typical 
problems cost only a few dollars to run. A reliable estimate of the accu- 
racy of the approximate $ is available. 

1.2 The Laplace package 

The method described in this paper has been implemented in the 
Laplace program package. A user's guide for the Laplace package, with 
several examples, is available separately. 10 The program package is 
written in EFL, 11 an extended Fortran language. The output of the EFL 
compiler is portable Fortran. 

The package solves Laplace's equation in two dimensions or in three 
dimensions with axial symmetry. Two-dimensional regions must be 
bounded by straight-line segments. Three-dimensional regions must 
be figures of revolution whose cross section in the (r,z) plane is bounded 
by straight-line, segments. The region may extend to infinity, but the 
boundary must not extend to infinity. The region may be multiply- 
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v 2 *(x,y) = fir + fir=o < pia ) 



connected. On each line segment, either $ or d$/dn may be specified 
as a boundary condition. In addition to smooth basis functions, the 
program automatically includes the appropriate singular basis functions 
greatly improving the achievable accuracy for regions with corners. 

Section II discusses the mathematical basis for the boundary integral 
equation method. Section III describes the implementation of the 
method. Possible extensions to the program package are discussed in 
Section IV. Section V has results for a sample problem. The appendix 
derives the new integral equation used. 

II. INTEGRAL EQUATION FORMULATION 

In this section, Laplace's equation is formulated as a pair of coupled 
integral equations. The two-dimensional partial differential equation 
is reduced to a pair of one-dimensional integral equations. 

We wish to solve 

d 2 4> d 2 <i> 
dx 2 dy : 

for (x,y) in a region D with boundary T. D may be multiply-connected, 
in which case T has several distinct parts. For now, we discuss only the 
two-dimensional "interior" problem, with D a finite region. At the 
conclusion of this section, we discuss the two-dimensional "exterior" 
problem, with D an infinite region, and the three-dimensional axisym- 
metric problem, both interior and exterior. 

As in Fig. 3, let (x s ,y s ) be the coordinates of the point at arc length s, 
and denote <£(s) = $(x«,y s ). We will use 4> for the potential of a general 
point, and <t> for a point on T. Let n s be the outward-pointing unit normal 
vector at s. For a point s not at a vertex of T, define 

\f/(s) = lim n s • V$(x,y). 

(x,y)—{x a ,y a ) 

The notation d<f>/dn s is also used for the right side of the above defini- 
tion. 

For the problem to be well-posed, a boundary condition must be given 
at each point of T. 12 The Laplace package allows the specification 

0( s ) = ^(s) on part of T, say 1^ (Plb) 

^( s ) = b 2 (s) on the remainder of I\ say T 2 . (Pic) 

For any fixed point F = (x F ,y F ), the Green's function, or fundamental 
solution to Laplace's equation, is 

G{x,y;x F ,y F ) = -V 2 ln [(x - x F ) 2 + (y - y F ) 2 ]. 

Except at point F, V 2 G{x,y;x F ,y F ) = 0. 
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Now let (x,y) be any point strictly inside D. Green's boundary identity 

3 

27r<l>(x,y) = P [iP(s)G(x s ,y s ;x,y) - 0(s)n s • V s G(x s ,y s ;jc,y)]ds. 

The gradient operator, V s , operates on the x s and y s . The above equation 
is usually abbreviated as 

2TT$(x,y) = £ h(s)G - <f>(s) |^1 ds, (1) 

with the arguments of G left implicit. 

If (x,y) is a point at arc length Una smooth part of T, it may be 
shown 3,13 that 

7r0(t) = ■£ I" *{s)G - 0(s) ^-1 ds. (2) 

The integral is now a Cauchy principal- value integral' at s = t. 

Suppose that the correct <f>(s) and \p(s) are not known, but only ap- 
proximate values 0*(s) and ^*(s) are known. Then the function $*(x,y) 
defined by 

27r<t>*(x,y) = f U*(s)G - 0*(s) — 1 ds 

exactly obeys Laplace's equation for (x,y) strictly inside D. <i>* will not 
obey the correct boundary conditions as {x,y) approaches the boundary 
unless (f)*{s) and ^*(s) are chosen correctly. 

Thus the boundary integral equation method is one of the class of 
"particular solution" methods. 14 - 15 Any approximate solution obeys the 
partial differential equation exactly, but only obeys the boundary con- 
ditions approximately. The advantage over the usual particular solution 
methods for Laplace's equation, as seen in Ref. 16, for example, is that 
the boundary integral particular solutions incorporate the exact 
boundary of the region and do not require a restricted region. They are 
more complicated to calculate, but are appropriate for the region. 

Equation (2) may be used to obtain an integral equation for 0* and 

r- 

dG l 
dn s ] 

Letting R be the vector from point t to point s, and R the length of R, 
(3a) may be written as 



tt0*(£) = -f \t*(s)G - 0*(s) ~ | ds. (3a l 



iii^^ * (s) ~ MW * (s) ] 



***(*) = + -^-0*(«) -ln(W*(s) ds (3b) 
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This boundary integral equation has been used for many years for solving 
Laplace's equation. 5-9 

For example, if </>(s) = b^s) is given on all of r (Dirichlet problem), 
set </>*(s) = 0(s), and the above is then an integral equation for the un- 
known \p*(s). Thus a two-dimensional partial differential equation has 
been reduced to a one-dimensional integral equation. An approximate 
solution may be obtained by expanding \f/*(s) in an appropriate set of 
basis functions, and taking a finite number of these. 

f («) = E ajfj(s). 

The /'s are piecewise constant functions in Ref . 6 and piecewise quadratic 
in Ref. 7. We discuss an appropriate set of /'s later. The integral equation 
(3a) then becomes 

£ a, -f fj(s)G ds = -f bM^-ds + xbiit). 
y=i Jr Jr drc s 

If M = N points ti are chosen at which to make this equation hold exactly 
(collocation), a set of N linear equations for the N unknowns, ay, is ob- 
tained. If M points ti, M > N, are chosen, an over-determined set of 
linear equations is obtained for the ay's. This reduces the sensitivity of 
the approximate solution to the exact choice of the £,. The equations 
are 

N 

£ A u aj = r b j = 1,2,. . . ,M, 



where 



A ij =-Cf j {s)G(s,t i )ds 

r i =-£b 1 (s)^^ 1 ds + T r b l (t i ). 

J r dn s 



These may be solved in a least-squares sense, say, by a standard sub- 
routine. 17 

In addition to the obvious advantages of this formulation, there is a 
well-known disadvantage. The integral equation for \p*(s) is a Fredholm 
integral equation of the first kind, 18 of the type 



u(x) = £ H(x,y)v(y)dy, 



where u and H are known and v is to be determined. This kind of integral 
equation is ill-conditioned (sometimes called ill-posed); it is difficult 
to obtain accurate solutions for j/. 19 > 20 The reason for the difficulty is easy 
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to see. Since v appears only inside the integral, its high-frequency com- 
ponents are not well-determined; by the Riemann-Lebesgue lemma, 



lim ( H(x,y) sin (ny)dy = 0, 

n— Jo 



if H is not too badly behaved. The difficulty numerically is that the 
matrix [Ajk \ is ill-conditioned. Small errors in calculating elements Ajk 
or rj lead to much-magnified errors in the coefficients a^. If a sequence 
of approximate solutions with increasing N is done, the larger matrices 
are increasingly ill-conditioned. The typical failure mode is that \{/*(s) 
does not converge as N increases, after a certain point; rather, spurious 
and unphysical oscillations in ^*(s) are seen. 

Various methods of ameliorating the difficulty have been suggested, 
such as regularization 21 and matrix singular-value decomposition. 19 ' 20 ' 22 
Better yet is to derive a Fredholm integral equation of the second 
kind. 

In the appendix, we derive the following identity for t any point at a 
smooth part of r. 

**(t) = -f U(s) f^ - [*(«) - *(t)] ~^- 



ds. (4a) 



For the Dirichlet problem, this identity leads to a Fredholm integral 
equation of the second kind for ^*(s) and is not ill-conditioned. Ap- 
parently, but surprisingly, (4a) is new. The integral equation derived 
from (4a) may be written as 

7r^*(t) = -^{[0*(s)-0*(O] 

v 2(n s • R)(n t ■ R) -fi 2 n s ■ n, n, ■ R i 

X — + R2 t*(s) ds. (4b) 

With two integral equations, one well-conditioned for 0* and the other 
for \p*, problem (PI) can be reduced to a set of linear equations with a 
well-conditioned matrix. If fitting point tj is on r x , where 0(s) is speci- 
fied, use (4). If tj is on r 2 , where \f/(s) is specified, use (3). A coupled pair 
of linear integral equations results. Analogously to the Dirichlet problem 
discussed earlier, appropriate basis functions and fitting points can be 
chosen, and the problem reduced to a set of linear equations. Some of 
the complications will be covered in later sections of the paper. 

For the Dirichlet problem, </> given everywhere on T, (4) cannot be used 
everywhere. Since (4) is independent of the zero of potential, (4) alone 
will lead to a singular matrix, of rank N — 1. Special methods may be 
used for dealing with rank-deficient matrices, or the other equation, (3), 
may be used at some fitting points. 
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2. 1 Exterior two-dimensional problems 

We now consider solving Laplace's equation in an infinite region, D, 
exterior to a finite boundary, T. To have a unique solution, it is insuffi- 
cient to specify either 0(s) or \p(s) at each point of T. In addition, the 
behavior of $(x,y) far from T must be specified. Use standard polar 
coordinates, (r,0), with r = Vx 2 + y 2 and suppose 

lim <*>(x,y) = ^ In - + *„ + 0(l/r), 
r —a» 2tt r 

where 4? m and ^«, are constants. If ^ m is specified, then a unique solution 
can be found. 12 ^a, is the negative of total flux extending to infinity. 

The earlier equations apply with small changes. For example, (3a) 
must be replaced by 

T**(t) = 2*$: + -f U*(s) | t*(s)G~\ ds, 

and the unknown <i>L must also be found. 

The user specifies ^ ro as well as boundary conditions. The Laplace 
package calculates an approximate value for $1 as well as for 0* and ip* 
onT. 

2.2 Three-dimensional axisymmetric problems 

Most of the preceding two-dimensional analysis needs only minor 
changes for the three-dimensional axisymmetric problem. Unlike the 
two-dimensional problem, the same formulation is adequate for interior 
and exterior three-dimensional problems. We use standard cylindrical 
coordinates (r,6,z). We wish to solve 

„,, „ . Id/ d4>\ 1 d 2 4? d 2 $ n 

in an axisymmetric region D; D is formed as a figure of revolution by 
rotating a region D, with boundary T, about the z-axis. The boundary 
conditions must also be independent of d x 

For any fixed point F = {xf,5>f,zf)> the Green's function is 

G(x,y,z;xF,yF,ZF) = J, , 9 , , x? 7 , T^TZ - 7T 

[(x - x F ) 2 + (y - y F ) 2 + (z - z F ) 2 \ l > 2 R 3 

The integral equation corresponding to (3) is 13 

27T0*(t) = -f- T^*(s)G - 0*(s) — 1 da 

The integral is an area integral on the boundary, the surface of revolu- 
tion. 
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Without loss of generality, let t be at 6 = 0. Let T be the intersection 
of fi with any plane 6 = constant. Express the area integral as an iterated 
integral over 6 and s, arc length along I\ Since <f>* and \J/* are independent 
of 6, all the 6 dependence is in G and dG/dn s . R 3 may be expressed as 

A|-Jti-2r.r*(l + coB9) l 

where 

With much manipulation, the integration may be performed, giving 
2*<t>*(t) = -C Ur a K(m)f*(s) 

+ [ »2 S " S ' R ' £(™) + 2n s • e r [K(m) - H;(m)]~U*(s)j £- , (5) 

where e r and e z are unit vectors in the = plane, and R2 is the vector 
in the = plane from t to s. 

R2 = (r s ~ fi)e r + (z s - z t )e z 

±r s r t 

rtm 

K and E are the usual complete elliptic integrals. 23 

E(m) = C * 2 (1 - m sin 2 u) l ' 2 du 

K(m) = C * 2 (1 - m sin 2 u)~ l ^du. 

The integral equation corresponding to (4) is considerably more com- 
plicated. 

r 4r f T E K — El 

2mP*(t) = -f--f n, • R 2 - 2r s n t • e r k*(s) 

JrRlll 1-m m J 

r £ 

+ (n s • e r )(n t • e r )(2E - K) - n s • n f 

L 1 — m 

(n s • R 2 )(n t ■ R 2 ) / 2(2 - m) v 

'" fl 2 (l-m) V(l-m) / 



_2 



+ — [(n s • e r )(n, • R 2 )r t - (n, • e r )(n s • R 2 )r s ] 
VI — m m I \ 



ds. (6) 



K and JS have been used as abbreviations for K(m) and E(m). 
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2.3 Error estimates 

If (x,y) is strictly inside D, the approximate potential (in two di- 
mensions) obeys 

2ir<f>*{x,y) = f \if>*{s)G - 0*(s) — 1 ds. 

If (x,y) is on a smooth part of I\ 

Tr^*(x,y) = -C- \yp*(s)G - **(*)— \ds- 

If (x,y) is at a vertex of I\ the ir is replaced by the interior angle of the 
vertex. Similarly, V<£* inside D and d"t>*7dn on T may be obtained by 
the analog of (4). 

The function $* as defined above exactly obeys Laplace's equation 
inside D. Therefore, by the maximum principle, 24 the maximum error 
in <£* occurs somewhere on T. 



$*(x,y) — $(x,y) ^ max 



**(x„y 8 ) - 4>(sY 



For the Dirichlet problem, a rigorous error bound is in principle possible 
by finding the largest discrepancy between <I>* and the boundary con- 
dition 0. However, finding the error bound can be more expensive than 
solving the integral equation. 

For mixed boundary conditions, a rigorous bound is in general im- 
possible. The above bound is still correct, but is not useful, since the true 
<i> is not known on all of the boundary. For certain restricted regions, 
another rigorous error bound can be obtained. 25 For these restricted 
regions, 



#*{x,y) - #foy) 



< max 
ri 



$*(x s ,y s ) - 4>(s) 



+ Rd max 

r 2 



f(s) 



where Rd is the maximum perpendicular distance from any point of 1^ 
to any other point of I\ This bound is in principle possible to compute, 
but is expensive in practice. 

An error estimate is available at no extra cost in the Laplace package, 
because of the method of solution. The over-determined system of linear 
equations is solved in a least-squares sense, minimizing the total fitting 
error (tfe), 



[M N "11/2 

E L (Ay-ay - r,) 2 , 
i=i y=i J 
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and returning this error. For the implementation discussed in the next 
section, with M « 3N/2, numerical experiments indicate that the TFE 
is a reliable upper bound on the error in the potential on the boundary, 
and a substantial overestimate of the error away from the boundary. 



III. IMPLEMENTATION 

Section II was general and discussed mathematics; we now become 
more specific and discuss numerical analysis. We also discuss some of 
the myriad details necessary to make a computer program feasible. 

3.1 Geometry 

Section II considered regions of arbitrary shape. The current imple- 
mentation of the Laplace package requires that T be composed of finite 
straight-line segments. This is in contrast to the usual practice in anal- 
ysis, of requiring that T be a smooth curve everywhere. Many practical 
problems have corners in their geometries, so it is essential to be able to 
handle such boundaries. It is easier to analyze exact corners than 
"smooth" geometries with a very small radius of curvature rather than 
a corner. In the following, each of the straight-line segments is called a 
side. 

3.2 Basis functions 

In the previous section, the choice of the basis functions fk (s) was left 
arbitrary. However, the particular choice made strongly affects the ac- 
curacy and efficiency of the program. At least four factors should be 
considered. 

(i) The basis functions should be able to model the behavior of <t>(s) 

and \J/(s) with only a few functions. 
(ii) If enough basis functions are used, they should be able to ap- 
proximate <t>(s) and \{/(s) arbitrarily well. 
(Hi) The basis functions should be a well-conditioned set, so that small 
errors in doing the integrals do not lead to large errors in the 
approximate solution. 
(iu) The integrals of the basis functions times G, dG/dn s , and 
d 2 G/dn s dn t must be tractable, either analytically or numeri- 
cally. 
Historically, (iu) has been dominant. Symm 6 approximated curved 
boundaries by straight-line segments and used piecewise constant basis 
functions. Hayes 7 allowed boundaries to be straight-line segments or 
arcs of circles, and used piecewise quadratic basis functions. Blue 9 al- 
lowed straight-line boundaries and allowed piecewise polynomial basis 
functions. For all these choices, the integrals in (3) and (4) can be done 
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(P. 7) 




7 = 



Fig. 4 — Local polar coordinates used for eq. (7). 

analytically. The above authors did not treat the axisymmetric problem, 
but the integrals in (5) and (6) are intractable analytically. 

If the boundary has corners, \f/(s) can be infinite at the corners, and 
piecewise polynomial basis functions will not be able to approximate \f/(s) 
well. For example, in Fig. 4, suppose the boundary conditions are as 
shown, and the interior angle is y m . Choose polar coordinates (p,y) 
centered at the vertex, with angle y = on the = side; let s = s at the 
corner. Then for small p, <f> has an expansion (in two dimensions) 26 

*(p,7) = L C n p"n sin (a n y), (7) 

n=l 

where 

(n - V 2 )7r 



« n = 



7. 



Thus $ is identically zero on the line 7 = 0. On the line 7 = y m , the 
normal derivative is identically zero. 

On the line 7 = 0, where s ^ s, the normal derivative is 

iA(s)= £ a n C n (s - s )»«-\ 

n=l 

For the case y m > 7r/2, we have ai < 1, and we expect i/'Cs) to be infinite 
at the corner, unless C\ happens to be zero. For \[/(s) to be approximated 
accurately with only a few basis functions, one of them should be «i(s 
— So)" 1-1 with the correct a\. 

Similarly, on the line 7 = y m , where s < so, the potential is 

0(s)= f C n (s -s)""(-l) n+1 . 

n=l 

Therefore a basis function (s — s)" 1 is needed on the 7 = y m side. In fact, 
only a single unknown coefficient, C\, need be introduced to deal with 
the worst part of 0(s) and \p(s) at s . It may also be desirable to include 
a few of the less singular basis functions. The Laplace package includes 
singular basis functions for which a < a max , with a default value a max 
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= 1, and does not include any singular function whose a is within 0.1 of 
an integer. 

Similar singular basis functions are used at corners where <p is specified 
on both sides and where tp is specified on both sides. Since the expansion 
(7) is not necessarily convergent far from the vertex, the singular basis 
functions centered at a corner are used only on the two sides meeting at 
that corner. 

For axisymmetric problems, the expansion (7) does not hold, but the 
exponent of the singularity is the same for off-axis points, since the 
singularity depends only on the highest order derivatives in the differ- 
ential equation. 26 For singularities on the axis, the exponents are dif- 
ferent, 27 and singular functions have not yet been implemented. 

Additional "smooth" basis functions are also needed. A well-condi- 
tioned family of basis functions is B-splines. A brief description of some 
of their properties follows. 28 - 29 B-splines are defined on a line divided 
into intervals by knots. B-splines of order k are piecewise polynomials 
of degree k — 1. Each B-spline is nonnegative, has exactly one maximum, 
and has local support. The sum of B-splines at any point is identically 
one. In any interval, exactly k B-splines are nonzero; each B-spline is 
nonzero in at most k intervals. The B-splines on a line are uniquely de- 
termined by the knots, which may be multiple. At a knot with multi- 
plicity m, a &th-order B-spline has k — m — 1 continuous derivatives. 
If m = k — 1, the B-spline is only continuous; if m = 1, the B-spline has 
k — 2 continuous derivatives. 

The user chooses a k, the same for all sides, and the number of interior 
knots on each side. Mesh spacing proceeds according to the following 
rules. A vertex is called singular if its expansion, as in (7), has «i < 0.9. 
If a singular basis function is used, the vertex is called compensated. If 
the two vertices delimiting a side are each either compensated or 
nonsingular, the interior knots on the side are spaced uniformly. Oth- 
erwise, the interior knots are spaced closer together near uncompensated 
singular vertices. 

For a given mesh, higher-order B-splines are potentially more accu- 
rate, since the approximation error can be 0(h k ) [30], where h is the 
maximum mesh. However, the integrals for higher order splines are more 
difficult, and there are more unknown spline coefficients for higher k. 
Currently, the Laplace package restricts k to be 2, 3, or 4. 

3.3 Boundary conditions 

If <f>(s) or \f/(s) is specified as an arbitrary function, the integrals in- 
volving the boundary conditions require special methods. Instead, the 
Laplace package does a least-squares fit of the boundary condition to 
a B-spline. A separate fit is done on each side; the same order and mesh 
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are used as specified by the user for that side for the unknown \p*(s) or 

0*(«). 

For the Neumann problem — \p specified on all of T — the problem is 
undetermined up to an additive constant in *. A solution exists for the 
interior problem only if f r tis) ds is exactly zero. The Laplace package 
currently will not solve the Neumann problem; must be specified on 
at least one side. 

3.4 Integrals 

With polynomial basis functions and boundaries composed of 
straight-line segments, the integrals in (3a) and (4a) can be done ana- 
lytically. This can cause conditioning problems; B-splines are a well- 
conditioned basis only if calculated properly. 29 The integrals in (5) and 
(6) cannot be done analytically, even with polynomial basis functions. 
With singular basis functions like (s - s ) a , none of the integrals can be 
done analytically unless a is special. 

All the necessary integrals can be done accurately and efficiently by 
the numerical methods described in this section. We first consider (3a) 
and (4a), for any fixed t. For straight-line boundaries, the integral over 
T is divided up into a sum of integrals over the line segments. We con- 
sider only a single segment, and eliminate any subscript referring to the 
segment. On the segment, n s is constant, and R may be written as 

B. = R ± n s + (s — s x )e 8 , 

where e s is a unit vector along the side. Also expand n t as 

n, = nj.n s + n\\e s . 



The portions of (3a) and (4a) from the segments are 



- V 2 ln [R\ + (s- s ± ) 2 ty*(s)\ds (3c) 

2R ± [n ± R ± + n\\(s - s ± )] - [Rj + (s - s ± ) 2 ]n ± 
[i?i + (s-s x ) 2 ] 2 



n ± R± + n\\(s — Sj_) 



R- ± + (s-s ± ) 2 



ds. (4c) 



We first consider the case where the point t is not on the segment in 
question. Then the Green's function parts of the integrals are not sin- 
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gular; either R ± ^ or s ± is not in the interval [si,s 2 ]. As an example, 
look at the integral with the logarithm in it, and look at one term of the 
expansion of ^*(s), with the basis function B,(s). Further divide the 
segment into subintervals, between knots of the spline, so that over each 
subinterval fij(s) is a polynomial. Over each subinterval, the integrand 
is a polynomial times a nonsingular function. Gauss-Legendre quadra- 
ture is ideal for such integrands, if the order of the quadrature rule can 
be determined a priori. (An automatic quadrature method could be used, 
such as Refs. 31 or 32, but these are usually less efficient.) Computing 
the order of the quadrature rule necessary can be done using the results 
of numerical experimentation. If s yi and s ;2 are the ends of the sub- 
interval, lets c = (sj 1 + Sj 2 )/2 and h = Sj 2 — s ;r The change of variable 
u = 2(s — s c )/h changes the term in question to 

- 7 f 1 [In (/i 2 /4) + In [a 2 + (u - 6) 2 ]]fi 4 (s c + hu/2)du, 
4 J-\ 

where a = 2R ± /h and b = 2(s ± — s c )/h. Since the integral is over a single 
mesh interval of B„ we expect the error to be no worse than the worst 
error in any of the fcth-order B-splines with fc-fold knots at —1 and 1, and 
no interior knots, since the latter B-splines vary more rapidly over the 
interval. Thus we look only at the errors in these k B-splines; call them 
B{u) to distinguish them. 

Now consider the family of integrals 

Ij(a,b,k) = f 1 In [a 2 + (u - b) 2 ]Bj(u)du. 

Let Ej(a,b,k,n) be the error in evaluating Ij(a,b,k) by an n point 
Gauss-Legendre quadrature rule, and let 



E(a,b,k,n) = \ £ Ej(a,b,k,n) 2 ] 



1/2 



Numerical experiments show that in the (a, 6) plane, the locus of constant 
E(a,b,k,n) is approximately an ellipse. For given n and desired accuracy, 
c, there is an ellipse with semi-axes A(k,n,e) and B(k,n,e) so that the error 
is satisfactory if a and b are outside the ellipse, or 



i— — T+r— — l 



2 



For doing the integrals Ij(a,b,k) to accuracy e, the functions A (fc,n,c) 
and B(k,n,e) are determined experimentally for a series of values of n. 
(The default values are n = 4, 6, 8, 10, 12, and 16, and e = 10" 6 .) For any 
particular a and b, the smallest satisfactory n is used. If the largest n 
available is insufficient, then the interval is divided; this is seldom 
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necessary. In the Laplace package, A{4,n,t) and B(4,n,t) are used for 
fc<4. 

The other integrals in (3a) and (4a) are done similarly, using numer- 
ically derived ellipses for B-spline basis functions. For singular basis 
functions, Gauss-Jacobi quadrature formulas are used; these are Gauss 
quadrature formulas on (0,1) with weight function x° _1 . Different 
quadrature formulas are used for each of the unique a's used in singular 
basis functions. The same ellipses as calculated for B-splines are used; 
slightly smaller ellipses could be used, but the gain in efficiency is 
small. 

The Gauss quadrature formulas are calculated portably by the method 
of Sack and Donovan, 33 using programs in the PORT library. 34 

The integrals of (5) and (6) are somewhat more complicated than those 
of (3a) and (4a), but are no harder numerically. The same ellipses are 
used. 

If point t is on the line segment, then the Green's functions in the in- 
tegrals have singularities. The integrals (3c) and (4c) simplify somewhat, 
since n s = n t , n|| = 0, n ± = 1, and R ± = 0. The 0* term in (3c) is iden- 
tically zero; the \p* term is 



-V 2 C 2 \n[(s-t)*]+*(s)ds. 

%Js\ 



The \t/* term in (4c) is identically zero; the </>* term is 



"f 



s Voo-0*(o] (/ " 



-*)»■ 



Special care must be taken to get accurate approximations to these 
singular integrals. 

First consider a B-spline basis function, B,(x), again dividing (si,s 2 ) 
into subintervals. If t is not in the subinterval in question, then the 
previous methods are adequate. (The Laplace package never takes t to 
be exactly at a knot.) For the logarithmic integral, the subinterval in- 
cluding t is, for some positive 5i and 5 2 , 

\n[(s-t)2]Bi(s)ds 

t-5i 

= - f ' In (u)Bi(t - u)du - f 2 In (u)J9,-(t + u)du 

Jo Jo 

= - 5 X f In (p)Bi(t - »Si)dv - In ($i) ]Bi(t-u)du 

- 8 2 f In {i>)Bi(t + vb 2 )dv - In (5 2 ) f ' Bi(t + u)du. 
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The integrals with In (v) are done by Gauss quadrature with weight 
function In (v); the others are done by Gauss-Legendre quadrature. The 
Gauss quadrature formulas with logarithmic weight function are also 
calculated portably by the method of Ref. 33. 

The Cauchy principal-value integral, for the subinterval including 
t,is 



"f 



»+*i B,-(s)-B,-(t) 

t-8! (S - 2 

«+*i ds r t+s 2 B i (s)-B i (t)-(s-t)B' i (t) 



Ut-h\ s — t Jt-bi 



ds. 



is - t)* 

The first integral is done analytically. The second has no singularity at 
s = t and is done analytically as 

This is adequate for low-order B-splines. For high-order B-splines, more 
care would be necessary. 

Now consider integrals with singular basis functions and with t on the 
same segment as s. If \f/(s) is given on the segment, 0*(s) may have (s — 
Si) a or (s2 — s) a terms; however, if ^(s) is given, (3) is always used for t 
on the side, and the <p* terms vanish because R j_ = 0. If 0(s) is given on 
the segment, \p*{s) may have (s — Si)" -1 or (S2 — s) a_1 terms; however, 
if 0(s) is given, (4) is almost always used, and the \J/* terms vanish because 
R j_ = 0. The exception, when 0(s) is given on a segment and (3) is used, 
occurs only for the Dirichlet problem, <j> given on all of T. Then (3) is used 
at the central fitting point of each side, and we need integrals of the 
form 

±V 2 a f 52 In [(s - t) 2 ](s - »!)«-l(fe. 

As much as possible of the integral 

>t 



s 



In {t - s)(s -si) a_1 ds, 



starting from s\, is done using Gauss- Jacobi quadratures. The remainder 
has only a logarithmic singularity. It and the integral form t to s 2 are done 
by Gauss quadrature with a logarithmic weight function, as described 
earlier in this section. 

3.5 Complete elliptic Integrals 

The complete elliptic integrals K(m) and E{m) are necessary for the 
axisymmetric problem. Suitable expansions are 35 
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K(m) = P K (1 -m)- Q k (l - m) In (1 - m) 

E(m) = P E (l - m) - Q E (1 - m) In (1 - m). 

Polynomial approximations for the Ps and Qs are given in Ref. 35. The 
argument (1 — m) is used instead of m to avoid excessive error asm—* 
1, i.e., as s —*■ t in (5) and (6). 

The combination D(m) = [K(m) — E{m)]/m is also needed. Asm-* 
0, D{m) —*■ 7r/4. For D{m), another approximation of the above type was 
generated. 



3.6 Fitting points 

At least N fitting points are needed to determine the N unknown 
coefficients of the basis functions. The work to calculate the matrix is 
proportional to the number of points used. If more than N points are 
used, the sensitivity of the solution to the placement of the points is di- 
minished, as is the amplification of any small errors in calculating matrix 
elements. In the Laplace package, approximately / times N fitting points 
are used; the default value of/ is 1.5. In the subinterval between each 
pair of knots, the number of fitting points is / times the number of un- 
knowns associated with the subinterval, rounded up. The fitting points 
are uniformly spaced within each subinterval. 

3.7 Scaling, constraints, and matrix solution 

Each row of the matrix corresponds to applying either (3) or (4) at one 
fitting point, £,-. To keep the solution approximately independent of the 
scaling of the region, each row corresponding to (4) is multiplied by the 
length of the side containing £,. 

For an interior problem, JV*(s) ds = 0. For an exterior two-dimen- 
sional problem, Ji/-*(s) ds = ^ m . Either restriction may be written as 
a linear equality constraint on the unknown coefficients. When the 
matrix equations are solved by QR factorization, such linear constraints 
can easily be enforced using a method described by Lawson and Han- 
son. 36 

3.8 Portability 

A portable stack allocation mechanism 34 is used for all temporary 
storage. The program is written in EFL. 11 The output of the EFL compiler 
is portable Fortran. 

Two parts of the program are not portable. For e ^ 10 -B , new ellipses 
are necessary. The approximations to the complete elliptic integrals are 
accurate to about 10 -8 . 
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IV. POSSIBLE EXTENSIONS 

In this section, we discuss several extensions to the Laplace package 
which could be implemented if there were sufficient incentive, and the 
difficulties involved with each. Combinations of the individual extensions 
pose further difficulties, but will not be discussed. 

4. 1 Higher-order B-splines and higher-accuracy Integrals 

B-splines of order higher than 4 are useful if very high accuracy so- 
lutions are desired. The only change necessary to allow higher-order 
B-splines or higher-accuracy integrals is to change the ellipse semi-axes. 
This feature was not included in the Laplace package, since calculating 
the ellipses portably for any specified accuracy and B-spline order re- 
quires too much code. Alternate methods for doing integrals to any 
specified accuracy are under consideration. 

4.2 Singular basis functions for axlsymmetrlc problems 

At a vertex away from the axis of revolution, an expansion similar to 
(7) will hold; the exponents [ot n \ are the same as for the two-dimensional 
problem with the same shape as the cross section of the figure of revo- 
lution. At a vertex on the axis of revolution, the exponents are different; 
on-axis singular functions have not been implemented. 

4.3 General linear boundary conditions 

In some applications, it is desirable to solve Laplace's equation with 
the general linear boundary conditions on T 

a(s)<fi(s) + b(s)+(s) = c(s), 

with a and b simultaneously nonzero. Then (3a), say, would become 

The difficulty here is in choosing a method for accurately evaluating the 
integrals involving a/b and c/b, unless alb and c/b are restricted dras- 
tically, say, to being constants on each of the boundary line segments. 

4.4 Curved boundaries 

Many applications have part or all of the boundary as a smooth curve, 
which the user might not wish to approximate by straight-line segments. 
In principle, all that is needed to allow r to be any smooth curve is a 
parameterization of x s , y S) and n s as a function of s. Again, the difficulty 
is in doing the integrals accurately and efficiently. The ellipse method 
would not be directly applicable. In addition, some of the integrals which 
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Fig. 5— A region which could be handled more easily by breaking it into two regions with 
an interface. 

vanished identically for straight-line boundaries would not vanish. For 
example, the <f>* term of (3b) for s and t on the same segment vanishes 
if the segment is a straight line, because n s -R is identically zero. If the 
line segment is curved, n s -R//? 2 in general has a finite limit as 8—t, and 
this term needs to be kept. 

4.5 Interlaces 

In other applications, there may be interfaces. A common problem 
is V 2 $! = in region D x with boundary I\, V 2 $ 2 = in region D 2 with 
boundary r 2 , and interface conditions on the common portions of Ti and 
IV Typical interface conditions are fa = fa and Kifa = K 2 fa, where *i and 
*2 are given constants. 

Implementing this extension would require a significant change in data 
structure, but otherwise would be easy. No new types of integrals would 
arise. The singular basis functions at corners which are also points on 
the common boundary depend on m and k 2 as well as the angles. 37 

This extension would also be useful for some single-region problems. 
A typical example is Laplace's equation inside a U-shaped region, Fig. 
5. This could be broken artificially into two regions as shown, with in- 
terface conditions <f>i = fa and fa = fa. The full region requires ap- 
proximately 3 N 2 integrals, if N is the number of unknowns. The two 
half-size regions would each have N/2 unknowns, plus a few extra for 
the boundary values on the dotted line. Each region would require 
somewhat more than 3(iV/2) 2 integrals, so that the total work would be 
somewhat more than half. 
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<l- =0 



•I' =0 



J/=0 



\fr=0 



* =1 tf= 

Fig. 6 — Geometry for a sample problem. 

V. AN EXAMPLE 

Figure 6 gives the geometry and the boundary conditions for a sample 
problem. The boundary conditions are \{/ = on the light sides, $ = 1 on 
the bottom dark side, and </> = on the top (L-shaped) dark side. This 
problem was solved approximately with third-order and fourth-order 
B-splines, and with various numbers of interior knots. Each side had the 
same number of interior knots. Some of the information is summarized 
in Table I. N is the number of basis functions. The running time is given 
in seconds for a Honeywell 6070 computer. TFE is the total fitting error, 
as defined in Section II. $\pds is over the side from 1 to 2. The next col- 
umn gives # at vertex 3. The approximate $ at (%,%) is the final column; 
vertex 1 is at (0,0) and vertex 3 is at (2,2). The time goes approximately 
as N 2 ; most of the work is in calculating elements of the matrix. Solving 
the matrix takes time proportional to N 3 , but the proportionality con- 
stant is smaller than that of the N 2 term. Figure 7 is a log-log plot of TFE 
against N; TFE seems to be converging as N -3 for third-order splines and 
as N -4 for fourth-order splines. These rates of convergence are the op- 
timum rates for approximating smooth functions by B-splines, 30 it is of 
interest to see them apparently applying for nonsmooth functions. The 



Table 



kord 


nknots 


N 


time 


TFE 


S+ds 


0at3 


*(%,%) 


3 





15 


1.19 


0.1067 


0.997242 


0.5482 


0.6189 


3 


1 


21 


1.98 


0.0255 


1.000287 


0.5069 


0.6196 


3 


2 


27 


3.21 


0.0231 


0.999932 


0.5011 


0.6193 


3 


3 


33 


4.91 


0.0123 


0.999951 


0.4999 


0.6192 


3 


4 


39 


7.09 


0.0080 


0.999965 


0.4997 


0.6192 


4 





21 


1.85 


0.0235 


0.999838 


0.5041 


0.6192 


4 


1 


27 


2.89 


0.0135 


0.999958 


0.4997 


0.6194 


4 


2 


33 


4.58 


0.0066 


0.999956 


0.4987 


0.6193 


4 


3 


39 


6.52 


0.0034 


0.999961 


0.4995 


0.6192 


4 


4 


45 


9.29 


0.0024 


0.999968 


0.4995 


0.6193 
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Fig. 7 — Total fitting error (TFE) vs the number of basis functions, for third -order and 
fourth-order splines. Lines proportional to N~ 3 and to N~ 4 are shown for comparison. 



last three columns appear to have converged to the accuracy allowed by 
the finite precision of the calculations. Matrix elements are calculated 
to a relative precision of about 10~ 6 , and the matrix has a condition 
number on the order of a few hundred, so accuracy of a few parts in 10 4 
is all that can be expected for boundary values. f\f/(s) ds can be more 
accurate, since the integration can average out the boundary errors. 

For these examples, the same number of knots was used on each side. 
Other examples may require differing numbers of knots on different 
sides. The intuition of the user is valuable in deciding on the number of 
knots per side. 



APPENDIX 



Derivation of Integral Equation (A) 

Equation (4) can be derived in various ways. We use a derivation 
modeled on the derivation of (3) as sketched in Ref. 13. Start with 
Green's identity in two dimensions. 



f f (uV 2 v-i>V 2 u)dA= C (u 



dv 



— v ) ds. 

drio/ 
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Fig. 8 — Integral equation (4). 



This identity is usually stated to hold for u and v which are C 2 inside D 
and C 1 in D + r, but is more generally true. For example, the condition 
on u can be weakened to include u's which have corner singularities as 
discussed in the body of the paper. The region D need not be simply- 
connected. 

Pick any fixed point (x t ,y>t) at arc length t on I\ at a smooth part of 
T. Let n t be the outward-pointing normal at t. Choose 

u(x,y) = *(x,y) - ${x t ,yt) 

v{x,y) = n t • VG{x,y;x t ,y t ), 

and apply Green's identity to the region D', which is D minus a sector 
of a circle, with radius e, centered at t (Fig. 8). Let T' be the circle sector. 
In D', V 2 u = and V 2 v = 0, so the area integral is zero. 

Consider the V integral, and use polar coordinates (r,6) centered at 
t. Let e r be the unit vector at (r,0) pointing away from the point t. On 
r, n s = - e r , 

v = —n t - e r /e, 
df/dris = — n, • e r A 2 . 
Expand <b(x,y) about (x t ,yt)- For (x,y) on T', 

$(x,y) = <S>(x,y t ) + ee r • V$(x t ,y t ) + 0(* 2 ), 

T^ (W) - ~*T- -^ ■ V*(x t ,y t ) + 0(e), 
dn dr 

where V<b(x t ,yt) is an abbreviation for V<f>(o;,y) | Xti y t . The integrals over 
T' may be evaluated explicitly in the limit as e-*-0. 
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lim Cu — ds = lim f* [ee r • V*(x t ,y t ) + 0( f 2 )][-n, • e r /e 2 ]edd 
t— o Jr dn s «— o Jo 

= - lim f * [e r • V*(x t ,y t ) + 0(e)][n t • e r ]dd 
t ^o Jo 

= - | n, • V*(x t> y t ) = - 1 M). 
In performing the integral, we used the identity 

J»7T 7J- 

(a • e r )(b • e r ) dd = — a • b, 
o 2 

true for constant vectors a and b. The other integral is evaluated simi- 
larly. 

lim I (— v ) ds 

e-o Jr \ dn s / 

= lim r \2tl*~\ [_e r . V^(x t ,y t ) + O(e))ed0 

= ->>. 

Thus the integral over T' gives —ir\p(t), in the limit e— *•(). Again in the 
limit €-*0, the integral over the remainder of T becomes a Cauchy 
principal-value integral, and (4) is obtained. The argument depends only 
on the most singular terms in the Green's function, and so is easily 
generalized to the axisymmetric case. 
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